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Polymeric thin films of various thicknesses, confined between two repul- 
sive walls, have been studied by molecular dynamics simulations. Using the 
anisotropy of the perpendicular, Pn(z), and parallel components, Pt(z), of 
the pressure tensor the surface tension of the system is calculated for a wide 
range of temperature and for various film thicknesses. Three methods of deter- 
mining the pressure tensor are compared: the method of Irving and Kirkwood 
(IK), an approximation thereof (IK1), and the method of Harasima (H). The 
IK- and the H-methods differ in the expression for Pt{z) {z denotes the dis- 
tance from the wall), but yield the same formula for the normal component 
Pn(z). When evaluated by MD (or MC)-simulations Pn(z) is constant, as 
required by mechanical stability. Contrary to that, the IKl-method leads to 
strong oscillations of Pn(z). However, all methods give the same expression 
for the total pressure when integrated over the whole system, and thus the 
same surface tension, whereas the so-called surface of tension, z s , depends on 
the applied method. The difference is small for the IK- and H-methods, while 
the IKl-method leads to values that are in conflict with the interpretation of 
z§ as the effective position of the interface. 
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The aim of statistical mechanics is to relate macroscopic quantities to microscopic degrees 
of freedom. An example for this connection is the virial equation of the pressure. Consider a 
system of volume V with M particles which interact by a pair potential U. Let the distance 
between two particles be denoted R (R = \R\). The pressure can then be written as a sum 
of two parts, 
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a kinetic (ideal-gas) part k-^Tp (p = M/V), which arises from the average kinetic energy 
and the momentum transfer of the particles on the container walls, and a potential part 
which accounts for the intermolecular forces. Two particles experience an interaction force 
— RU'(R)/R. When weighing the corresponding virial, —RU'(R), with the average density, 
p( 2 \R), of a particle at distance R from another one and integrating over all possible separa- 
tions, one obtains the contribution of the potential to the pressure. There are different ways 
to derive Eq. ([[J) (see [0-0], for instance), but none of these routes can readily be generalized 
to inhomogeneous systems. They all use the isotropy of space somewhere in the derivation 
and take p as a scalar. 

In inhomogeneous systems, however, the pressure in general depends on the spatial di- 
rection and on the position r where it is determined: It is a tensor P(f). Nonetheless, the 
pressure tensor can still be split into a kinetic part, P^, and a potential part, P u : 

P(r) = P K (r) + P u (r) . (2) 

The kinetic part may be expressed by a generalization of the ideal-gas contribution, 

P K (r) = k B Tp(r)l } (3) 

where p(r) is the density at r and 1 a 3 x 3 unit matrix. 

On the other hand, there seems to be no unique expression for P u (r) |2],|^-[10| . The 



origin of this problem may be explained as follows: The pressure tensor can be defined by 
the infinitesimal force dF acting across an infinitesimal surface dA which is located at r: 

dF(r) = -dA ■ P(r) . (4) 

If a particle moves across dA, the resulting momentum transfer contributes to P K (r). Since 
the momentum is associated with the particle position, it is a single particle property which 
may be well localized in space (see however [0). The ambiguity in the calculation of P(r) 
arises from the interaction between two particles: Which particles should contribute to the 
force at r? Somehow the non-local two-particle force, —U'(R), has to be reduced to a local 
force d-F(r) f7j. This ambiguity was already pointed out in the seminal work of Irving and 
Kirkwood, and they required that "all definitions must have this in common - that the stress 
between a pair of molecules be concentrated near the line of centers. When averaging over a 
domain large compared with the range of intermolecular force, these differences are washed 
out, and the ambiguity remaining in the macroscopic stress tensor ([0) is of negligible 
order" (footnote on p. 829 of @). 

In the present paper, we apply common ways to calculate P u (r) to a model of a glassy 
polymer film and determine the surface tension as a function of temperature. This work 
serves as a preparation for simulations on the sluggish relaxation of the film in the super- 
cooled state |i~3| . The paper is organized as follows: In Sect. |I| we discuss the theoretical 
background of various approaches to P u (r). Section [TTI| presents details of the model and 
simulation technique, and Sect. |TV| compiles the results. The final section contains our 
conclusions. 



II. THEORETICAL BACKGROUND 
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A. The Methods of Irving and Kirkwood and of Harasima 

Irving and Kirkwood || gave a definition of the P^-tensor by starting from a statistical 
mechanical derivation of the equations of hydrodynamics and by making a special choice for 
the particles that contribute to the local force: Only those pairs of particles should give rise 
to dF(r) for which the line connecting their centers of mass passes through the infinitesimal 
surface dA (see Fig. [l]) ||. With this choice they obtained the following expression for the 
potential part of the pressure tensor 

p^(r) = -i J ™ U\R) (j da p (2 \r - aR,r + (1 - a)R)j d 3 R , (5) 

where RR is a dyadic, U'(R) = dll/dR, and p( 2 \r; r') denotes the two-particle density 

pV\r-r , ) = (^5(ri-r)5(r J -r')) . (6) 

Using Eq. (^) one obtains from (g) 

pC7 ( r ) = -^(E ! ^ f/ 'M fdaS^-r + ar^)), (7) 

where = Vj - 7*j (r„ = |ry|). 

Equation ([5]) can be interpreted as follows: The term —RRU'/R is a tensorial general- 
ization of the virial —RU' of the integrand in Eq. ([!]). It accounts for the force RU' /R that 
a particle at r\ experiences from another particle at r 2 (R = r<i — r{). The virial has to 
be multiplied by the probability of finding two particles at r\ and r-i- The probability is 
proportional to the density p^tVi; r 2 ) which depends explicitly on both particle positions 
for inhomogeneous systems. Therefore, different values of p^ 2 \ri] r 2 ) are obtained for fixed 
R when shifting particle 1 or 2 to position r, where the pressure shall be determined, i.e., 
for t*i = r (a = 0) or r 2 = r (a = 1) (see Fig. 0). The integral over a takes all of these 
contributions into account. The outer integral finally sums over the possible vectors R which 
pass through dA. Equations @ and (0) are general and apply to systems of any shape if 
the particles interact by a pair potential. In the following we are interested in thin (poly- 
mer) films confined between two impenetrable walls. For systems with planar geometry the 
pressure tensor, P, depends only on the distance, z, from the wall |p|,[10|. Furthermore, the 
non-diagonal components of P vanish in thermal equilibrium and it can be written as (see 
Sect. [TDD 

P(z) = (e x e x + e y e y )P T (z) + e z e z P N (z) , (8) 

where e x , e y , e z are orthogonal unit vectors and the lateral, Pt(z), and normal component, 
Pn(z), of P(z) are given by 

PZZ{Z)=M*) alld P XX (Z) = Pyy(z) = P T ( Z ) . (9) 

Using 
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1 ~ I ' Z — Zi\ ^ / z 



da 5{z — azij — zi) = -, — r6 ( ) 6 



I Z{j | V Zij / V Zij 



and averaging Eq. ([?[) over the tangential coordinates one obtains fTOjJIi 



P u ( z ) = I / / P u ( r )dxdy 



A 



-^(E^^'M^e^e^)), (io) 



Tij | 2^" | \ Zij / V 



where A is the area of a plane in tangential direction. With Eq. (|K]) this leads to the 
following (full) expressions for the normal and tangential components of the pressure tensor 
for planar systems (IK-method) 

»-^^ W (^)e(^)) , (ii) 
*M*-i(E^e(^)e(^)) , (12) 

where p(z) denotes the density at z averaged over tangential coordinates x and y. These 
equations are valid only in thermal equilibrium (for an extension to non-equilibrium situa- 
tions see HH). 

In addition to the IK-expressions the formulas of Harasima are often used in the literature 
They are obtained from a different choice of the contributing interactions (see Fig. |T]): 
Harasima considered a prisma whose base is dA. The force dF(r) is thought to result from 
all interactions between particles in the prisma and those on the side of dA to which the 
vector dA points. This also includes particles whose center line does not pass through dA. 
Harasima's choice corresponds to a contour which goes parallel to the walls (or the planar 
surface) from r\ to (x2, V2,z\) arid then along the normal to T2 0,0- Using this convention 
he obtained the same results for the normal component as Irving and Kirkwood [Eq. ([H])], 

P*(z) = Ft*(z) , (13) 

but a different expression for the lateral component of the pressure tensor |||5| 

P*(z) = p(z)k B T --^(]T i-^rV,,) 5(z t -z)). (14) 

Thus, the tangential component, Pt, of the pressure tensor is not uniquely defined. Con- 
sequently, the pressure anisotropy, P^ — Pt, is ambiguous. This ambiguity is extensively 
discussed in the literature P,|4|- |lT)| . [T4| 



However, the integral over z of Eq. ([121) is identical to that of Eq. (14). This implies that 
both the IK and the H-methods yield the same results for any physical quantity which does 
not depend on the local profile of the pressure tensor. In particular, they lead to the same 
values of the surface tension 7 (Kirkwood-Buff formula M) 
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P n (z)-Pt(z) 



dz 
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r 2 - ^7 2 
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(15) 
(16) 



The factor 2 arises from the existence of two walls at z = —D/2 and z = D/2 in our simu- 
lation, D being the distance from one wall to the other (i.e., the film thickness). However, 
moments of Pn — Pt, such as the so-called "surface of tension" z s , i.e., the position where 
the surface tension acts, 



1 

27 



-D/2 



D/2 



Pk(z)-P t (z) 



dz 



(17) 



depends on the different choices made to determine P u . 
Harasima ||. 



This was already pointed out by 



In Sect. |TV| we want to show for the polymer model considered that the differences in 
z s obtained from the IK and H-expressions are small compared to the size a of a particle, 
but not negligible. The ambiguous nature of z s was discussed in detail in [pUT0|. In |I| a 



liquid-vapor interface is studied. Since there are no density oscillations near a free surface, 
which are characteristic of liquid- wall interfaces |P],|15|, we expect the difference between the 
IK and H-expressions for P-r(z) to be more pronounced for the thin films studied here. 



B. The Method of Planes 

Todd, Evans and Daivis |5]|5| have introduced a variant of the original IK-derivation to 
determine the pressure tensor (termed "method of planes") which avoids the ambiguity of 
defining a contour to relate two interacting particles. The problem is, however, not circum- 
vented because one has to choose a gauge for both the pressure tensor and the momentum 
density M. The derivation starts from the continuity equations for the mass and momentum 
and leads to 



M 

~2A\ 

i=l 

i+3 



for the potential part of the pressure tensor and to 



^(*) = t(E^^-^)) ( 20 ) 

i=l 



for the kinetic part (a = x, y, z), where m is the mass of a particle. In Eq. (18) sgn(x) is 
the sign function (= 1 if x > and —1 for x < 0), and F a i is the a-component of the force 
exerted on particle % by all other particles. Furthermore, 0(x) denotes the Heaviside step 
function and p a i is the a-component of the momentum of particle i. Using the identity 



5 



Q(zi — z) Q(z — Zj) — Q(zj — z) Q(z — Zi) 

one can verify that the diagonal components of the Eqs. ([H]) and (|20D yield the IK-expression 
for the normal pressure [Eq. (P]) j. Since Eq. (|T8|) contains a single sum instead of the double 
sum of Eq. (]TT|) , it is computationally more convenient. Therefore, we used Eqs. (|18|) and ( |20|) 
to calculate the normal pressure. However, these equations are not sufficient for determining 
the surface tension 7, as they do not contain the diagonal components of the pressure tensor 
parallel to the walls, i.e., P xx and P yy . On the other hand, they provide a method for the 
calculation of the viscosity ||. 



C. An approximate formula: IKl-method 



In the literature (see [|16|,|17]], for instance) there is still another formula for the pressure 
tensor, which is a kind of a "tensorized" version of the Harasima expression flliD (called 
"IK1" in U) 

P IK1 (.) = p(z)k B T i - ±( V r Jpi U'( rij ) 6( Zi -z)). (21) 



r 



'j 



Todd, Evans and Daivis || noticed that Eq. (^T|) is equivalent to a zeroth-order approxima- 
tion of the (full) IK-expression and gave a physical interpretation of the approximation in 
k— space (see Eq. (24) in ||]). One can also find a real-space interpretation in the following 
way. If one replaces the integral over a in Eq. ([7]) by the value of the integrand at the lower 
bound a = 0, one obtains 

=-\(Y, rj P l U '^) ^ ~ r )) ' ( 22 ) 



which gives the potential part of the IKl-expression (^1|) after averaging over the tangential 
coordinates. 

Thus, the IKl-method corresponds to the assumption that the two-particle density 
p^(r 1 ;r 2 ) is unchanged upon translation of both arguments along the line R = r 2 — r\ 
which connects the points 1 and 2. However, the breaking of translational invariance is one 
of the basic characteristics of inhomogeneous systems. The more the system is inhomoge- 
neous, the more the IKl-expression ( pl|) for P^(z) should become inaccurate. On the other 
hand, integration over z yields the same result as the IK- and H-approaches. Therefore, the 
IKl-method leads to the same surface tension 7, but to a different value for z s compared to 
the other two methods. 



In Sect. [TV] we show that the IKl-result for z s is too large to allow for an interpretation 
of z s as the effective position of the interface, i.e., as the distance of closest approach of a 
particle to the wall. Furthermore, Eq. ( pl~D leads to strong oscillations of Pn i n contrast to the 



condition of mechanical stability which requires a constant profile for Pn (see section [11 D|). 
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D. Mechanical Stability Requires Pn = const 



In equilibrium, mechanical stability requires that the gradient of the pressure tensor 
vanishes 

V ■ P = , (23) 

where denotes the null vector. For a system with planar symmetry, the non-diagonal 
components of P must also vanish (otherwise shear forces would exist) and the lateral 
components should be identical. So, we have 

dP XX dP m OP; 



dx Bx + ~df ey + ~dz~ Bz = and p ^ = p yy^- ( 24 ) 

Since dP xx /dx = 0, dP yy /dy = on the one hand, and P xx = P yy on the other, the lateral 
components can be functions of z only. Furthermore, since dP zz /dz = 0, the normal 
component of the pressure tensor is independent of the distance from the surfaces and must 
be identical to the external pressure PN.ext- This gives 

Ps(z) = P zz = ^N.ext = const and P T (z) = P xx {z) = P yy (z) , (25) 

i.e., Eq. (H). The argument presented is not new. It essentially follows the discussion of |2| 
(see p. 44 of [Q]). We repeated it here to stress the erroneous character of expression (pip . 
In Sect. we will see that only the IK- (or H-) formula flllD satisfies condition (p5|) . The 
independence of Eq. ( O ) on z was already proved analytically in the work of Harasima (see 
p. 224 of H). This important properties helps us to set the pressure in the simulations for 
a given wall separation and temperature. 



III. SIMULATION OF POLYMERIC FILMS 



A. Model 



We study a Lennard- Jones model for a polymer melt |TJ| embedded between two impen- 
etrable walls. All simulation results are given in Lennard- Jones (LJ) units. Two potentials 
are used for the interaction between particles. The first one is a truncated and shifted LJ- 
potential which acts between all pair of particles regardless of whether they are connected 
or not, 



tWr) 







-U LJ (r c 



if r < r c , 
otherwise , 



where 



U LJ (r) = 4e (r/a) 12 - (r/a) 



and r c = 2 x 2 1 / 6 . The connectivity between adjacent monomers of a chain is ensured by a 
FENE-Potential g9| 
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^feneO") = ---Ro ln 



1 - 



r 

R 



2i 



where k = 30 is the strength factor and Rq 
The wall potential was chosen as 

U w (z) = 



1.5 the maximum allowed length of a bond. 



a 



(26) 



where z = | ^particle — ^waiil (^waii = ±-D/2). This corresponds to an infinitely thick wall 
made of inifinitely small particles which interact with inner particles via the potential 
180(r/cr)~ 12 /(7rp wa n) where p wa n denotes the density of wall particles. The sum over the 
wall particles then yields (a/z) 9 . 

The static and dynamic properties of this model were studied in the bulk when gradually 
supercooling towards the glass transition [18,20-2^]. The model begins to develop sluggish 
relaxation if the temperature drops below T rs 0.7 and yields a critical temperature of mode- 
coupling theory of T c bulk ~ 0.45 |20| upon further cooling. We quote this value for the sake 
of comparison with the film results to be discussed below. 



B. Contribution of the Walls to the Normal Pressure 

As the wall potential acts only in normal direction, the expressions (|T2| ) and ([14]) for 
remain unchanged. To obtain the contribution of the walls to Pn one can consider each wall 
as an additional particle of infinite mass and use Eq. ([18]) for the extended system of M + 2 
particles. Starting from Eq. ( ]T5| ) one can show that 

M 

^botwall , 



pwallsJK^ = ^ ^ F W (Zj - Zbotwall) ©0*i ~ z) Q(z 
i=l 
1 M 

_ T\ ^w(^topwall — z i) ©(^topwall — z ) ®(z ~ z i)j > (27) 



where Fw(z) = — dUw(z)/dz, ^botwaii < Z{ < -2 top waii for all (inner) particles (i.e., excluding 
the wall particles) and Zbotwaii < z < z t o P waii for all planes. From Eq. (|27D it follows that the 
force Fw of a wall on a particle contributes to the normal pressure on a given plane if the 
plane lies between the particle and the wall. 

Similarly, one can derive the contribution of the walls within the IKl-approximation by 
starting from Eq. fl2~l"|). This yields [p4| 



pwalls,IKl^ = /y^ Fw ( z . _ 2botwaU ) S(zi - Z) 
1=1 
1 M 

- -A ^ ^w(^topwall - Zi) S(Zi - Z)^ , (28) 
i=l 

where the sum runs over inner particles only, as before. Since F w (zi— z')5(zi — z) is equivalent 
to Fw(z — z')5(zi — z), P^ alls ' IK1 (z) can be written as a product of the density profile and a 
contribution from the walls, i.e., 



S 



P™ ' 1 (z) = F W (Z- Zbotwall) - -Fw(^topwall ~ z) p(z 



C. About the Simulation 

The equilibration of the system was done in the NpT-Ensemble. The production runs, 
however, were performed in the NVT-Ensemble because we are also interested in analyzing 
the dynamics of the films later on (for preliminary results see [(R|). 

At the beginning of the simulation the velocities of all particles were set to zero and 
NRRW- (Non- Reversal- Random- Walk-) chains were "synthesized", i.e., only the average 
bond length and bond angle (known from previous bulk simulations) were used to build a 
chain of N (= 10) monomers. This initial state corresponds to very high energies (usually 
E(t = 0) > 10 10 ) due to the occurrence of extremely short distances between non-bonded 
monomers. 

The surplus of energy must be removed to prevent numerical instabilities. For the bulk 
this can be done by replacing the full LJ-potential by a softer one. The LJ-potential is then 



switched on smoothly [19]. For our model, however, it was necessary to keep the (full) wall 
potential from the very beginning of the simulation to avoid penetration of the walls. We 
thus left the potentials unchanged, but used an adaptive time step: First, the maximum 
force -F max and the maximum velocity v max were determined. A time step A was then chosen 
so that the resulting displacement of a particle, which is subject to -F max and moves with 
initial velocity t> max in direction of -F max , would be dr max = 10~ 3 . This (empirical) value is 
only applicable if -F max does not point in direction of a bond vector whose size b is closer 
to the maximum bond length Rq (see ?7fene) than 10~ 3 , since a displacement of this size 
could break the bond. In such a situation we chose dr max = (R — 6 max )/2 instead of 1CT 3 
to adjust the time step (6 max denotes the largest measured bond length). The equations of 
motion were then integrated with this time step and the procedure was repeated. 

After about 250 MD steps the velocities of all particles were renewed by drawing them 
from the Maxwell distribution, and the time derivative of the volume was set to zero. These 
steps are important to warrant the numerical stability of our procedure. Our criterion for 
the end of this stage was that the minimum distance between particles should not be smaller 
than a certain value, empirically 0.8, and that the normal pressure of the system should not 
be too far away from the external value, i.e., |-Pn(0 — PN,ext\/PN,ext < 10~ 2 , where P\j(t) 
was computed as an average over the last 20 samples preceding time t. The sample distance 
was empirically chosen to 10exp(l/T) MD steps to take into account stronger correlations 
at lower temperatures. Since we kept the film thickness D fixed, the simulation at constant 
pressure was realized by varying the area (= A) of the simulation box parallel to the walls. 
During this initial stage a high bath temperature, T = l, was used. 

After this initial stage (with a typical duration of 10 5 MD steps) the time step could 



be set to A = 0.003. This value is close to that used in previous bulk simulations JTS 
The system was then slowly cooled down to the desired temperature by gradually reducing 
temperature in a step-wise fashion: The bath temperature was set to the next smaller value 
and the system was propagated for a a certain amount of time before the bath temperature 
was decreased again. 
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At the end of the cooling process the sampling of the mean-square displacement of the 
chain centers parallel to the walls, g3\\(t), and of the volume was started. The system was 
propagated until g> 3 || > §R 2 ee \\ where P ee || denotes the component of the chain's end-to-end 
vector parallel to the walls. This criterion suffices to reach the free diffusive limit and to 
equilibrate the system completely. During this period the system volume was sampled once 
every 1000 time steps and the average volume of the system was calculated. The equilibrated 
configuration was then further propagated until the instantaneous volume reached the aver- 
age value within a given relative accuracy, usually 10~ 5 . At this point the program fixes this 
volume and switches to a (pure) Nose-Hoover-Algorithm (NVT-Ensemble) for production 
runs in the canonical ensemble. During a production run sampling was done once every 
1000 time steps. 



IV. RESULTS 

A. Profiles of P N (z): IK1 versus (full) IK 

In order to analyze the pressure profiles for our model we studied different film thicknesses 
(P = 3,5, 10,20) at various temperatures while always keeping PN,ext = 1- For this external 
pressure many results for the bulk behavior are known |IEpU| -|2"3|l . Here, we want to discuss 
two representative cases: D = 3 (~ 2R g where R g ~ 1.45 is the bulk radius of gyration) 
at T = 1, and D = 10 (~ 7P g ) at T = 0.42. The temperature T = 1 corresponds to the 
high-temperature, (ordinary) liquid state of the melt, whereas T = 0.42 belongs to the 
supercooled temperature regime close to the critical temperature of mode-coupling theory 
(T C (P = 10) « 0.39 0). 



For a film of thickness D = 3 ten independent runs of 10 6 time steps were simulated 
at T = 1 and PN,ext = 1- The total number of particles was 1000 corresponding to 100 
chains of length N = 10 (this number of monomers per chain was always kept fixed in our 
simulations). For D = 10 five independent runs were done at T = 0.42. The length of a run 
was 4.4 x 10 7 time steps. Samples were taken every 1000 steps. The much longer simulation 
time in this case is necessary to allow for a detailed analysis of the dynamics of the system 
which is very slow at this temperature. 

Figures ^| and || show the simulation results for the normal component of the pressure 
tensor, Pn, calculated according to the IK- and IKl-prescriptions, respectively [see Eq. ( |Tl"D 
and Eq. ([211) 1 . Furthermore, they resolve the different contributions stemming from the 
kinetic part, the virial (forces between inner particles, i.e. excluding the walls) and the 
walls. The striking difference between both prescriptions is that the IKl-method yields 
strong oscillations, whereas the pressure profile of the IK-method is constant throughout 



the film, in agreement with the condition of mechanical stability [see Sect. [PDj. This 
deficiency of the IKl-method has already been pointed in an analysis of the pressure tensor 
for a simple liquid |8]||. 

Since the kinetic contribution to Pn is proportional to the density profile p(z), Fig. || 
shows that practically no particle is present in the vicinity of the walls. The excluded-volume 
interaction creates a depletion zone of about 0.8 between the wall (z wa ii = ±1.5) and the 
monomer positions at this temperature. Any plane in this region separates all particles of 
the system, which lie on the side of the plane facing towards the inner part of the film, from 
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the wall on the other side. There is no interparticle force across the plane and thus the virial 
contribution to the normal pressure vanishes. The behavior of P^(z) near the wall arises 
only from the wall-particle interaction. This interaction does not depend on the position of 
the plane as long as all the particles stay on the opposite side, i.e., as long as p(z) ~0. This 
explains why Pn is constant in the region close to the walls. With increasing distance from 
the wall the density starts to increase from zero. Then, the kinetic and virial parts begin 
to contribute, whereas the effect of the walls decreases. In this intermediate region none of 
the contributions is negligible, but their sum still remains constant, in accord with Eq. (p5|). 
Very far from the walls the contribution of the walls to Pn becomes negligible. There, one 
expects that the variations of the kinetic and virial terms must be opposite to each other. A 
first indication of this opposite behavior can be observed in Fig. |2|. A better demonstration 
is, however, shown in Fig. f| where the film thickness is large enough to exhibit an inner 
region with negligible wall contribution. 

Contrary to that, the various contributions of the IKl-methods are (almost) in phase. 
Figure [| illustrates that the strong deviation of P^ K1 from a constant is caused by the 
interaction of the wall with the monomers close to the maximum of p(z) if D = 3. If the 
film thickness increases, Fig. |5| shows that the oscillations of Pn propagate through the 
whole film. Close to the wall, the dominant contribution still comes from the wall-monomer 
interaction, whereas the oscillations in the inner part of the film are in phase with the 
virial. The contribution of the virial is negative close to the wall, reflecting a predominantly 
attractive interaction between the monomers. This dominance of the attractive interaction 
is also visible for the (correct) IK-method, but is much less pronounced in this case. 

The situation becomes more complicated when studying the lateral component of the 
pressure tensor. Here, the two alternative formulas, Eqs. (12) and (14), can yield completely 
different profiles. Figures |6| and ^ compare the IK and the H- versions to calculate the lateral 
pressure Pt(z) for D — 3, T— 1 and D— 10, T = 0.42, respectively. Whereas both methods 
oscillate in phase with one another for the thicker film, they are anti-correlated for D = 3. 
The lateral pressure of the IK-method is positive close to the walls, but negative in the 
middle of the film, whereas the behavior is just vice versa for the H-method. Due to the 
aforementionend ambiguity of Pt{z) it is impossible to decide which methods yields the 
physically more realistic result. If the film thickness increases, the qualitative difference 
between the IK- and H-methods (almost) vanishes and only quantitative differences remain. 
The oscillations of Pt(z) clearly reflect the monomer profile. In the inner portion of the film 
they are much weaker for the H-method than for of the IK-method. This is related to the 
local nature of Eq. (|T4T ) due to the presence of delta-function. Density oscillations are thus 
incorporated not only in the kinetic term, but also also in the virial part of the Harasima 
formula. Both terms partially cancel each other. Although the profile generated by Eq. (|14"D 
is thus closer to PN, ex t than that of Eq. fljjD , this should not be considered as an argument 
in favor of the H-method. A clear distinction between both methods would only be possible 
if one could find a quantity which specifically probes Pt(z) and whose behavior is known a 
priori, as it was the case for P^(z). 



11 



B. Surface Tension and Surface of Tension 



As mentioned in Sect. [II A| , integration of the pressure profiles over z yields the same 
result for the IK-, H- and IKl-expressions. Therefore, all methods must lead to the same 
surface tension 7 [i.e., Eq. (0)]. This expectation is nicely borne out by the simulation 
data for all film thicknesses and temperatures studied, where 7 was calculated by Eq. (p~5|) . 
Figure ||] exemplifies this behavior for D = 5 (~ 3i? g ). With decreasing temperature the 
surface tension increases by about a factor of 1.5. 

Qualitatively, this temperature dependence is expected. The monomer density of a 
polymer melt close to a hard wall exhibits a profile that is large at the wall and decays 
towards the bulk value in an oscillatory fashion with increasing distance from the wall (see 
Fig. 11 as an example) fl5| . Since the average density grows with decreasing temperature 
in a simulation at constant pressure, the maxima and minima of the profile become more 
pronounced. This means that there are more monomers in the highly populated layers at 
low than at high temperatures, and that the oscillations of profile become more long-ranged. 
These effects tighten the film so that the free energy needed to move monomers out of the 
interface, i.e., the surface tension, should increase as temperature decreases. The same effect 
is expected when reducing the film thickness because the layering is more pronounced in 
thinner films. This expectation is borne out by the simulation data (see Fig. |^). 

Contrary to 7, the discussion of Sect. [11 A| implies that the surface of tension, z B , depends 



on the method applied. This fact is illustrated in Fig. [10| which shows the temperature 
dependence of z s for the IK-, H- and IKl-methods. The difference between IK and the 
H-methods is rather small, whereas the IKl-result lies substantially above the values of the 
other two methods. Since z s can be interpreted as the distance of closest approach of a 
monomer to the wall, i.e., as the effective position of the wall, the following simple argument 
rules out the IKl-result: At temperature T, a particle can only penetrate into a (soft) wall 
up to the point, z w , where the wall potential balances thermal energy of the particle, i.e., 
U w {z w )/T=l. Using Eq. © this gives 

1 \ 1/9 

(29) 



T 

Equation ( p9|) is compatible with the IK- and H-predictions, but not with the IKl-result. 
Another way to illustrate this point is shown in Fig. |TT] where we plotted the monomer density 
profile of a film of thickness D = 10 at T = 0.42. With increasing film thickness the IK- and 
H-values for z w approach one another - for D = 20 for example, they are indistinguishable 
within the error bars (not shown here) -, but the disparity to the IKl-result remains. The 
figure clearly shows that the IKl-method places the effective wall position deeply into the 
interior of the film, whereas it has to be situated in the region where the density profile 
approaches zero. 



V. CONCLUSIONS 



We have reported simulation results for the pressure tensor of polymeric thin films which 
investigate the ambiguity in the definition of the potential part of this quantity. We studied 
three common methods: the method of Irving and Kirkwood (3| , that of Harasima f| and 
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an approximation of the IK- method, the so-called IKl-approach ||. On a microscopic scale, 
our simulation results show significant differences between the IK- and H-methods for the 
lateral component Pt(z) of the pressure tensor. However, both methods agree with each 
other for the normal component P^(z). They lead to a constant profile in accord with 
mechanical stability. On the other hand, the IKl-formula exhibits strong oscillations of 
Pn(z), as also found in [|,@. The origin of this discrepancy comes from the fact that IK1- 
method corresponds to a zeroth-order approximation of the IK-expression, which assumes 
translational invariance of the two-particle density p( 2 )(ri; r 2 ) with respect to the difference 
vector R = r 2 — V\. This assumption is not valid in thin films which exhibit density 
oscillations that are damped out only gradually with increasing distance from the wall. This 
local structure becomes more pronounced with decreasing temperature and film thickness. 
The more pronounced it is, the stronger the IKl-method will deviate from the IK-expression. 

However, when integrated over the whole system all methods give the same result. Thus, 
the surface tension, 7, of a planar system can still be calculated using each of these methods. 
This is no longer possible for moments of the pressure profiles, such as the surface of tension 
z s . The fact that IKl-expression can be used to calculate the surface tension although it 
is based on an incorrect expression for the local pressure tensor has occasionally caused 
confusion in the literature. For instance, Pandey et al. |T7[ applied the IKl-expressions 
to polymer films confined between one repulsive and one attractive wall, taking the local 
pressure profiles literally. The present analysis shows that the pressure profiles published 
in |17j are incorrect. Thus we hope that the present analysis will help to avoid this confusion 



in future simulation studies. 
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FIGURES 




FIG. 1. Schematic illustration of the different contributions to P (r) which are taken into 
account by Irving and Kirkwood (IK-method) and by Harasima (H-method). Let dA be an in- 
finitesimal surface situated at position r [panel (a)]. In the IK-method all particles whose center 
line passes through dA contribute to the force felt across the surface [panel (b)], whereas Harasima 
assumes that the interaction between the particles inside a prisma with base dA and those on 
the side to which dA is pointing causes the force at r [panel (c)]. Panel (b) shows two possible 
contributions in the IK-method. If R = ri — r±, the position vectors of the particles can also be 
expressed as: r\ = r — nR and r2 = r + (1 — a)R (0 < a < 1) [see Eq. (||)]. The interaction 
between r' x and r' 2 is also taken into account in the H-method, but not that between r\ and r<i- 
On the other hand, particles at r% and (= r% + R) contribute in Harasima's approach, whereas 
they don't in the IK-method. 
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FIG. 2. Different contributions to the normal pressure profile Pn(z) for a film of thickness 
D = 3 2i? g ) at T = 1 (high-temperature liquid state) and -PN,ext = 1 according to the (full) 
IK-method [see Eq. (0)]. The H-method yields the same result [see Eq. (|l3|)l. The various parts, 
kinetic (full line), virial (dashed line) and wall (dash-dotted), mutually balance one another to yield 
a constant profile Pn(z) = -PN,ext (circles), as required by mechanical stability (see Sect. II Dj ). The 
difference between PN,ext = l (vertical dashed line) and Pn(z) shows the accuracy to which we can 
fix Pn ext in the simulation for this film thickness. The difference is smaller than 2%. 
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FIG. 3. Different contributions to the normal pressure profile Pn(z) for a film of thickness 
D = 3 (w i? g ) at T= 1 (high-temperature liquid state) and -PN,ext = 1 (vertical dashed line) according 
to the IKl-method [see Eq. (H)]. Contrary to Fig. |, the various parts, kinetic (full line), virial 
(dashed line) and wall (dash-dotted), do not balance, but amplify one another, resulting in a 
(non-physical) oscillatory structure of Pn(z) (circles). 
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FIG. 4. Different contributions to the normal pressure profile Ptss(z) for a film of thickness 
D = 10 (w 7.Rg) at T = 0.42 (supercooled state close to T c w 0.39 p3|) and -PN.ext = 1 (vertical 
dashed line) according to the IK-method [see Eq. (|TT|)1. The H-method gives the same result 
[see Eq. (|i~3|)]. As in Fig. ^, the various parts, kinetic (full line), virial (dashed line) and wall 
(dash-dotted), mutually balance one another and sum up to a constant profile Pn(z) = -Pi\r,ext 
(circles), in agreement with the condition of mechanical stability (see Sect. |IID ). 
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FIG. 5. Different contributions to the normal pressure profile Ptss(z) for a film of thickness 
D = 10 (w 7R g ) at T = 0.42 (supercooled state close to T c « 0.39 |L3]]) and P^^xt = 1 (vertical 
dashed line) according to the IKl-method [see Eq. (^j])]. As in Fig. ^, the various parts, kinetic 
(full line), virial (dashed line) and wall (dash-dotted), give rise to a non-constant pressure profile 
(circles) contrary to the requirement of mechanical stability. 
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FIG. 6. Tangential component Pr(z) of the pressure tensor as obtained from the IK-formula 
[Eq. (|l2|)] and from the H-formula [Eq. fli~4|)] for D = 3 (~ 2i? g ), T=l (high-temperature liquid 
state) and PN,ext = l- The thin solid line shows the kinetic contribution k^Tp^z) (divided by 15 to 
put it on the scale of the figure). 
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FIG. 7. Tangential component Pt{z) of the pressure tensor as obtained from the IK-formula 
[Eq. (0)] and from the H-formula [Eq. (|l§] for D = 10 (rj 7-R g ), T = 0.42 (supercooled state 
close to T c « 0.39 and -PN.cxt = 1 (vertical dashed line). The thin solid line shows the kinetic 
contribution ksTp(z). 
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FIG. 8. Temperature dependence of the surface tension, 7, calculated by Eq. fllq), using the 
IK-, H- and IKl-methods for D = 5 (~ 3R g ) and PN,cxt = l- The temperatures shown range from 
the high-temperature, liquid state of the film to the supercooled state. 
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FIG. 9. Comparison of the temperature dependence of the surface tension, 7, for D = 5 
(~ 3i2g) and D = 20 (« 14R g ). The results of the IK-method are shown only. The other methods 
(H- and IKl-methods) yield the same 7's within the error bars, as exemplified in Fig. ^ for D = 5. 
The external pressure is Pjss ext = 1- The temperatures shown range from the high-temperature, 
liquid state of the film to the supercooled state. 
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FIG. 10. Temperature dependence of the surface of tension z s [Eq. (|i~7|)1 determined by the 
IK-, H- and IKl-methods for D = 5 and -PN,ext = 1- The solid line shows the simple estimate, 
Zw = I/T 1 / 9 [Eq. (||)], for the position of the wall. 
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FIG. 11. Monomer density profile of a film of thickness D = 10 (« 7R g ) at T = 0.42 
(> T c w 0.39 |l~3[| ) and i^N.ext = 1- Since the profile is symmetric around the middle of the 
film, the figure only shows one half of it. The scale of the abcissa was shifted so that the wall is 
placed at z = 0. The vertical lines mark the values of z s computed according to the IK-, H- and 
IKl-methods. 
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